//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Programa de Tomografia Roche para VCs e pVCs. 
// 
// - Versao Nov - 2008: 
//     Novas funcionalidades:
//       - Considera o fluxo de cada linha em separado.
//       - Pode excluir velocidade orbital no espectro modelo. Para aplicar em
//        dados corrigidos.
// - Versao out - 2009:
//    Corrigido falha na leitura do parametro de obscurecimento de borda
//
// Inclui biblioteca dataFileReader.h (Zabot A.) para ler arquivos de dados.
//
// Tiago Ribeiro - LinuxC++ 2008.
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#include <iostream>
#include <cmath>
#include <fstream>
#include <cmath>
#include <vector>
#include <sstream>
#include <string>
#include <cstring>

using namespace std;

#include "convolve.hpp"
#include "gsl_multimin_fdf.hpp"
#include "light_curve.h"
//#include "light_curve.cpp"
//#include "roche.cpp"
#include "../dataFileReader_v0.1.1/dataFileReader.h"
#include "rtomF.hpp"

//#include "lib_zabot/dataFileReader.h"


int main(int argc, char** argv)
{
  test_args(argc,argv);

	
  cout << "# - Loading data and preparing internal structures...";
  cout.flush();
  
	
  DataKeeper Rdata(argc,argv);
  Rdata.fchoice(Fslow_3);
//  Rdata.fchoice(Fdefault);
		

  //setando imagem inicial (I = 1.0) e parametros para inicializar a classe gsl_multimin_fdf
  gsl_vector *img0 = gsl_vector_alloc(Rdata.img_sz());


  for(int i = 0; i < Rdata.img_sz(); i++) gsl_vector_set(img0, i, 1.0);
  Rdata.mval(*img0);

	double fwhm = Rdata.FWHM();
	Rdata.Recalc_Default_Map(*img0,fwhm);

  double mfac = Rdata.mfac();
	double pfac = Rdata.S(*img0) > 0 ? -mfac*Rdata.X2_slow_3(*img0)/Rdata.S(*img0) : mfac*Rdata.X2_slow_3(*img0);

	//pfac = 1./pfac;
  Rdata.X2_reset();
	
  //pfac*=0.0005*pfac;
  
  Rdata.pfac(pfac);
 //setando parametros da gsl 
  const double STP_SIZE=1.0*sqrt(Rdata.img_sz());
  const double TOL = 1e-8;
  const double cfac1=2.0;
  gsl_multimin_function_fdf MinFunc;
  MinFunc.f = &objective;  // definicao da funcao
  MinFunc.df = &objective_div; // derivada
  MinFunc.fdf = &ob_ob_div; // calcula f e df junto!
  MinFunc.n = Rdata.img_sz(); // dimensao da funcao
  MinFunc.params = &Rdata; // Parametros
  
  cout << "[DONE]\n"
	<< "# - Initializing GSL_multimin_fdf Algorithm ...";
  cout.flush();
	
  GSL_Multimin_fdf MinObj(Rdata.img_sz(),MinFunc,img0,STP_SIZE,TOL);
  
  cout << "[DONE]\n";
  //int pfac_cont=0;
	
  //cout << Rdata.S(*img0) << "\n";

  double s2 = Rdata.Sc(); //(Rdata.Sc()+Rdata.S())/2.0;
  double x2 = Rdata.X2c();//(Rdata.X2()+Rdata.X2c())/2.0;
  double score2 = (x2-Rdata.Xaim()+1.)/(Rdata.Xaim()+1.);
	
	
	//ofstream path("s_x2_path.dat",ios_base::app);

	const int NITER = Rdata.nIter();
	
	// pre processing iteration
	/*
	cout << "# - Pre proc iteration ...";
	cout.flush();
	MinObj.iterate();
	cout << endl;
	pfac = Rdata.S(*img0) > 0 ? -mfac*Rdata.X2_slow(*img0)/Rdata.S(*img0) : mfac*Rdata.X2_slow(*img0);
	Rdata.pfac(pfac);
	MinObj.restart();
	 */
	cout << "# - Starting Penalty function: " << Rdata.pfac() << " / " << Rdata.Xaim() << endl << "# - Minimizing ...\n";


  while(MinObj.iteration() < NITER){  
	
    int status = 0;
	double s = Rdata.Sc(); //(Rdata.Sc()+Rdata.S())/2.0;
	double x = Rdata.X2c(); //(Rdata.X2()+Rdata.X2c())/2.0;
	double score0 =  MinObj.score();

	if( x <= Rdata.Xaim() )  {
	  cout << "# - " << Rdata.X2() << " " << Rdata.Xaim() << " " << Rdata.pfac() << endl
		<< "# - break (Xaim Reached)\n";
	  break;
	  }

	//while( ( x > Rdata.Xaim()+1 ) && ( MinObj.iteration() < NITER ) ) {// ?? oq colocar aqui?
	  while( ( x > Rdata.Xaim()) && ( MinObj.iteration() < NITER ) ) {// ?? oq colocar aqui?
    //for(int i = 0; i < 10; i++){
		cout << "## - Itr[" << MinObj.iteration() << "] started ";
		cout.flush();
      status = MinObj.iterate();
		cout << " [DONE] \n";
//	  double dif_S =  (s==0 ? 0. : fabs((s-Rdata.S(*MinObj.BestEstimate()))/s));
	  double dif_S =  (s==0 ? fabs(Rdata.S_gauss_fast(*MinObj.BestEstimate())) : fabs((s-Rdata.S_gauss_fast(*MinObj.BestEstimate()))/s));
		
		//if (s != 0.0) dif_S = fabs((s-Rdata.Sc())/s);
		cout << "## S[best] = " << Rdata.S() << " - S[curr] = " << Rdata.Sc() << endl
			<< "## X2[best] = " << Rdata.X2() << " - X2[curr] = " << Rdata.X2c()
			<< " | Score = " << MinObj.score()<< endl;
      cout << "##  DS = " << dif_S << " - DX2 = " << fabs((x-Rdata.X2c())/x) << " - DScore = " << fabs(score0 - MinObj.score())/score0 <<endl;
/*		cout << "## - S = " << Rdata.S(*MinObj.BestEstimate()) << endl
		<< "## - S = " << Rdata.S_gauss_fast(*MinObj.BestEstimate()) << endl
		<< "## - X2 = " << Rdata.X2_slow_2(*MinObj.BestEstimate()) << endl;
*/

		double bfac = 1e-2;
      //if( ((dif_S < bfac) || (fabs((x-Rdata.X2c())/x) < bfac)) && (fabs(score0 - MinObj.score())/score0 < bfac) )
	   // */
		
		
		if( ((dif_S < bfac) && (fabs((x-Rdata.X2_slow_3(*MinObj.BestEstimate()))/x) < bfac)) && (fabs(score0 - MinObj.score())/score0 < bfac) ){
			cout << "# - Storing values...\n";
			s = Rdata.S_gauss_fast(*MinObj.BestEstimate()); //Rdata.Sc(); //(Rdata.Sc()+Rdata.S())/2.0;
			x = Rdata.X2_slow_3(*MinObj.BestEstimate()); //(Rdata.X2()+Rdata.X2c())/2.0;
			score0 =  MinObj.score();		
			

			break;
		}

		
		
		if (dif_S > bfac) {
			cout << "# - Recalculating Default Map...";
			cout.flush();
			Rdata.Recalc_Default_Map(*MinObj.BestEstimate(),fwhm);
			//Rdata.Recalc_MLFLUX();
			cout << endl;
		}
		

		cout << "# - Storing values...";
		cout.flush();
		s = Rdata.S_gauss_fast(*MinObj.BestEstimate()); //Rdata.Sc(); //(Rdata.Sc()+Rdata.S())/2.0;
		x = Rdata.X2_slow_3(*MinObj.BestEstimate()); //(Rdata.X2()+Rdata.X2c())/2.0;
		score0 =  MinObj.score();		
		
		cout << "\n#\n";
		


		//if( Rdata.X2() <= 1.0){
		//	cout << "# - Rdata.X2 break \n"; 
		//	break;
		//}
	
    }
	
	  
	  cout << "###############\n## - s_old = " << s2 << endl 
		<< "## - s_new = " << s << endl 
		<< "## - Ds = " << ( (s2 > 0) ? (s-s2)/s2 : s ) << endl
		<< "## - x2_old = " << x2 << endl
		<< "## - x2_new = " << x << endl
		<< "## - Dx2 = " << (x-x2)/x2 << endl
		<< "###############\n";

	  if( ( ( fabs(( (s2 > 0) ? (s-s2)/s2 : s )) < 1e-6 ) || ( fabs((x-x2)/x2) < 1e-4 )  ) ){ //|| ( fabs(score2 - MinObj.score())/score2 < 1e-5 ) ) {
		  cout << "# -(s == s2) && (x = x2) break\n "; 
		  break;
	  }
	  
	  
	  s2 = s;
	  x2 = x;
	  score2 = MinObj.score();
	  
	  //cout << Rdata.Xaim() << endl;
    //break;
    double cfac0=cfac1; 
    //double cfac2=1.001;
    
    Rdata.pfac(Rdata.pfac()*cfac0);
	//Rdata.Recalc_MLFLUX();
	Rdata.Recalc_Default_Map(*MinObj.BestEstimate(),fwhm);
	cout << "# - New Penalty function: " << Rdata.pfac() << endl << "#\n";
    MinObj.restart();

	
  }

	ostringstream path_name;
	
	path_name << argv[1];
	string path_filename = path_name.str();
	
	string::size_type ph_frst = path_filename.find_first_of(".");
	string path_root (path_filename, 0, ph_frst);
	
	ostringstream ph_saida;
	ph_saida << path_root << "_path.dat";
	
	ofstream path(ph_saida.str().c_str());

	path << MinObj.iteration() << " " << Rdata.S_gauss_fast(*MinObj.BestEstimate()) << " " << Rdata.S(*MinObj.BestEstimate()) << " " << Rdata.X2c() << " " << Rdata.pfac() << endl;
	path.close();

  cout << "# - Minimization finished with " << MinObj.iteration() << " major iterations." << endl
       << "#\tX2 = " << Rdata.X2_slow_3(*MinObj.BestEstimate()) << endl
      << "#\tS  = " << Rdata.S(*MinObj.BestEstimate()) << endl
	<< "#\t S_g = " << Rdata.S_gauss_fast(*MinObj.BestEstimate()) << endl;
  ostringstream saida0;
  
  saida0 << argv[1];
  string filename=saida0.str();
  
  string::size_type frst = filename.find_first_of(".");
  string out_root (filename, 0, frst);
  
  ostringstream saida;
  saida << out_root << ".rimg";
  ofstream oo(saida.str().c_str());
  
  if(!oo) {
    ostringstream erro;
    erro << "Arquivo de saida nao pode ser aberto\n\tArquivo: " 
	<< saida;
    aborta(erro.str());
  }
  
  cout <<  "# - Writing resulting image in " << saida.str() << " ..." << endl;
  
  Rdata.show_par(oo);
  oo << "# - S_g = " << Rdata.S_gauss_fast(*MinObj.BestEstimate()) <<
	" | S = " << Rdata.S(*MinObj.BestEstimate()) << endl;
  
  surface ImgEnd = Rdata.image();

  Rdata.Recalc_Default_Map(*MinObj.BestEstimate(),fwhm);
	gsl_vector * dft_map = MinObj.BestEstimate(); //Rdata.Default_Map();
	
  for(int i = 0; i < Rdata.img_sz(); i++){
    
//        cout << "aqui 1\n";

    //double I = gsl_vector_get(MinObj.BestEstimate(),i); //img_res.at(i)/double(ga.populationSize()); 
    double I = gsl_vector_get(dft_map,i); //img_res.at(i)/double(ga.populationSize()); 
  //      cout << "aqui 2\n";
    //if(I < 1) I = 1;

    oo <<ImgEnd.theta[i] << " "
      << ImgEnd.phi[i] << " " 
      << ImgEnd.px[i] << " " 
      << ImgEnd.py[i] << " "
      << ImgEnd.pz[i] << " "
      << ImgEnd.nx[i] << " "
      << ImgEnd.ny[i] << " "
      << ImgEnd.nz[i] << " "
      << I << " " 
      << ImgEnd.vx[i] << " "
      << ImgEnd.vy[i] << " "
      << ImgEnd.vz[i] << endl;
  }

  oo.close();
  string ss = saida.str();

	
  Rdata.write_res(ss);
	
	
  ostringstream saida2;
  saida2 << out_root << "_dft.rimg";
  ofstream oo2(saida2.str().c_str());
	
	if(!oo2) {
		ostringstream erro;
		erro << "Arquivo de saida nao pode ser aberto\n\tArquivo: " 
		<< saida2;
		aborta(erro.str());
	}
	
	cout <<  "# - Writing resulting default image in " << saida2.str() << " ..." << endl;
	
	Rdata.show_par(oo2);
	
	
	ImgEnd = Rdata.image();
	
	Rdata.Recalc_Default_Map(*MinObj.BestEstimate(),fwhm);
	
	dft_map = Rdata.Default_Map();
	
	for(int i = 0; i < Rdata.img_sz(); i++){
		
		//        cout << "aqui 1\n";
		
		//double I = gsl_vector_get(MinObj.BestEstimate(),i); //img_res.at(i)/double(ga.populationSize()); 
		double I = gsl_vector_get(dft_map,i); //img_res.at(i)/double(ga.populationSize()); 
		//      cout << "aqui 2\n";
		//if(I < 1) I = 1;
		
		oo2 <<ImgEnd.theta[i] << " "
		<< ImgEnd.phi[i] << " " 
		<< ImgEnd.px[i] << " " 
		<< ImgEnd.py[i] << " "
		<< ImgEnd.pz[i] << " "
		<< ImgEnd.nx[i] << " "
		<< ImgEnd.ny[i] << " "
		<< ImgEnd.nz[i] << " "
		<< I << " " 
		<< ImgEnd.vx[i] << " "
		<< ImgEnd.vy[i] << " "
		<< ImgEnd.vz[i] << endl;
	}
	oo2.close();
	
  return 0;
}

